program bench_compute_y

  use compute, only : compute_y

  implicit none

  integer, parameter :: nloops = 1000
  integer, parameter :: nmeasuremax = 1488
  integer, parameter :: nmeasure = 1286

  integer :: i, yi, yj, count_0,count_rate, count_max
  real :: t0, t1

  real, dimension(nmeasuremax) :: stepsize_y, y_current, n0, reject_y_vector,dumx,dumy
  real, dimension(nmeasure)    :: z, m1,dum2
  integer, dimension(nmeasure) :: timeindex_nonzero
  real, dimension(nmeasure,nmeasure) :: Dinv
  real, dimension(nmeasuremax,nmeasuremax) :: Rinv_current
  integer :: reject_y
  logical :: in_subsample
  real    :: m0T,n0T

  write(*,*) "Allocating data..."

  reject_y = 0
  call random_number(m0T)
  call random_number(n0T)
  call random_number(stepsize_y)
  call random_number(y_current)
  call random_number(n0)
  call random_number(reject_y_vector)
  call random_number(dumx)
  call random_number(dumy)
  call random_number(z)
  call random_number(m1)
  call random_number(dum2)
  call random_number(Dinv)
  call random_number(Rinv_current)
  do i = 1, nmeasure
    timeindex_nonzero(i) = i
  end do
  
  write(*,*) "Starting benchmark..."
  call system_clock(count_0,count_rate,count_max)
  t0=real(count_0)/real(count_rate)

  ! assuming timeindex_nonzero is monotonic increasing
  yj = 1
  in_subsample = .false.
  do i = 1, nloops
     do yi=1,nmeasuremax
        if (yi .eq. timeindex_nonzero(yj)) in_subsample = .true. 
        call compute_y(yi,y_current,stepsize_y(yi),m1,z,nmeasure,&
                       dum2,Dinv,n0,dumx,dumy,Rinv_current,n0T,m0T,reject_y,&
                       reject_y_vector,timeindex_nonzero(yj),in_subsample,.true.)
        if (in_subsample) then
          yj = yj + 1
          in_subsample = .false.
          if (yj .gt. nmeasure) yj = 1
        endif
     enddo

  end do

  call system_clock(count_0,count_rate,count_max)
  t1=real(count_0)/real(count_rate)

  write(*,*) "Benchmark completed in ",t1-t0, "seconds"

  
end program
